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Conduction between graphene layers is suppressed by momentum conservation whenever the layer 
stacking has a rotation. Here we show that phonon scattering plays a crucial role in facilitating 
interlayer conduction. The resulting dependence on orientation is radically different than previously 
expected, and far more favorable for device applications. At low temperatures, we predict diode- 
like current-voltage characteristics due to a phonon bottleneck. Simple scaling relationships give a 
good description of the conductance as a function of temperature, doping, rotation angle, and bias 
voltage, reflecting the dominant role of the interlayer beating phonon mode. 



Graphene has generated broad excitement both for 
fundamental science and for its potential applications in 
technology Attention has turned increasingly to bi- 
layer and few layer graphene , because of their scientific 
richness and their promise in technological applications 
involving band-gap opening by an external electric field, 
high current-carrying capacity, and electro-optical cou- 
pling Interlayer conductance is crucial in almost 
all such applications. However, most fabrication methods 
lead to "twisted" layer stacking, with a random angle of 
rotation between layers. As a result, momentum conser- 
vation in certain respects "decouples" the layers fl3l — TjTqI ] . 
Strictly speaking, the layers do couple, as evidenced by 
velocity renormalization [2()l - [24l | . But the states near the 
Dirac point in one layer largely decouple from those near 
the rotated Dirac point in the other layer (l7Hl"9l |25| . 
suppressing interlayer current [i~9j . 

Understanding the interlayer conduction has proven 
surprisingly subtle and difficult. Calculations of the in- 
terlayer conductance to date have required some phe- 
nomenological lifetime broadening, with the conductance 
depending on this broadening [l9j . In the limit of small 
carrier scattering the incoherent transport picture will 
fail, and it becomes problematic to even define an inter- 
layer conductance independent of the external contacts 
[19 ] . With a realistic scattering lifetime the conductance 
is well defined, but it's calculated value is extraordinarily 
sensitive to even small details of the matrix-element mod- 
eling [lil, HH . Finite conductance is predicted even for 
twisted stacking [lj| due to Umklapp electronic coupling 
fl~8j , but the coupling deca ys e xponentially with the size 
of the rotational supercell 17 - HJ 2f| . (Here "supercell" 
denotes the primitive cell of the Bravais lattice of the 
commensurate rotated bilayer.) This would pose a se- 
rious obstacle to many device applications, where some 
reliable lower bound on the interlayer conductivity is es- 
sential. Experimentally, some interlayer transport is ob- 
served [261.1271] . but the mechanism remains unclear. Ex- 
trinsic scattering mechanisms could also relax momen- 
tum conservation, especially for defect-rich substrates. 

Here we show that essentially all of the conceptual and 
computational problems of twisted bilayers are resolved 
by including phonon scattering explicitly in the calcula- 



tion of interlayer transport. Then incoherent transport 
is assured, no phenomenological broadening is required, 
and the calculated conductance is not so sensitive to 
the matrix-element model. The phonon-mcdiatcd con- 
ductance arises primarily from the vibrational "beating" 
mode of the bilayer, in which out-of-plane acoustic (flex- 
ural) modes of the two layers are 180° out of phase, so 
the layers vibrate against each other. We introduce a 
simple tight-binding Hamiltonian that can describe the 
electron scattering by this mode for arbitrary rotation. 

The resulting dependence on rotation angle is quite 
different than previously expected, and more favorable 
for device applications. In particular, the extraordinary 
sensitivity to rotation angle disappears, replaced by a 
smooth and mild dependence; and the interlayer conduc- 
tance is never very small at room temperature. Thus 
device behavior can be robust despite random rotation 
angles. Simple scaling relationships give a good descrip- 
tion of the conductance as a function of temperature, 
doping, rotation angle, and bias. Also, at low tempera- 
ture we find a diode-like turn-on of the conductance with 
increasing voltage. 

The phonon mediated current between the layers is 
given by 28| : 
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where the sum is over k,k',s,s',fx and spin degener- 
acy has been included. Here k and s label 2D elec- 
tron wavevector and band, respectively, (i labels phonon 
bands, / is the Fermi function, and V is the voltage bias 
between the layers. To first order, 
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where n 



is the Bose-Einstein occupancy, and 
ksl-n'e-ph^k's') is the electron-phonon matrix el- 
ement for carrier scattering from state ^> ks to &k's' by 
emitting or absorbing a phonon hu qfi in branch [i with 
wavevector q = k — k' . The electron energy E k in each 
monolayer is described by the tight binding model with 
the in-plane first neighbor hopping 0] t = 3.1 eV; and 
only out-of-plane phonons are included in evaluating 
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i/e-p/i- The conductance is G = I/(VA), where A is the 
area. 
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FIG. 1. (Color online) Dependence of bilayer coupling and 
conductance on rotation. Calculated band splitting A is 
shown vs (a) rotation angle 0, and (b) the corresponding pe- 
riod L of the rotational supercell. Electronic conductance (i.e. 
calculated omitting phonon scattering) at T = 0, with doping 
level Ef — 0.26 eV in both layers j2j|], is shown vs (c) angle 
and (d) period L. See text for reference conductance Gj 
and for the solid line in (d). Phonon mediated conductance 
at room temperature and same doping is shown vs (e) and 
(f) L. The solid black curve in (e) is given by Eq. ()S1|) with 
only one adjustable parameter, g=0.34 eV/A. Deviations at 
angles < 2kp/K ~ 3° are expected, see text. 



To render the calculation computationally tractable 
requires a simple tight-binding model such as that of 
Ref. [lo| . That model or any reasonable model would suf- 
fice for demonstrating the qualitative effects of clcctron- 
phonon coupling, which can be directly understood from 
the requirements of momentum and energy conservation 
in the electron+phonon system. However, to obtain rea- 
sonable quantitative accuracy, we desire a model that 
correctly describes the dependence on both rotation and 
interlayer spacing, insofar as these are known from LDA 
calculations. Several authors have examined this pro blem 
under the constraint of fixed interlayer spacing [2lL [2fj| , 
which is sufficient in the absence of phonons. The "band 
splitting" A in the bilayer [l7| is often used to conve- 
niently represent the degree of interlayer Dirac^point cou- 
pling, and hence the interlayer conductance [3, . (A 
is the energy from the second valence band at K to the 
second conduction band, shown in Fig. 3 of Ref. [l8|.) We 



require that the model also reproduces the few available 
LDA results [24[ for the dependence of A on layer spac- 
ing. This rules out the existing models, with the possible 
exception of that in Ref. (25|, which is too cumbersome 
for our purposes. We therefore adapt the simplest avail- 
able model [30( with a phenomenological lateral screening 
length \ xy , 
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where h a = 3.35 A is the equilibrium interlayer distance, 
Tij is distance between carbon atoms in adjacent layers, 
and £ij = {{xi — Xj) 2 + (j/j — Uj) 2 ) 1 ^ 2 is the in-plane dis- 
tance between the atoms. Here we take t± = 0.4 eV. We 
obtain a reasonable fit to the available LDA results for 
A in rotated bilayers [I?} and dependence on interlayer 
spacing 24j using \ z = 0.6 A, \ xy = 1.7 A and a=1.65. 
The model also reasonably well reproduces LDA results 
for the Fermi velocity renormalization (24[. 

The phonon eigenvalues and eigenvectors are calcu- 
lated using an atomistic valence force model (3lj plus 
a Lennard- Jones (LJ) 6-12 potential [32|, 33 1 between 
atoms in adjacent layers. The clcctron-phonon Hamil- 
tonian is obtained by expanding Eq. <|3j> to linear order 
in the phonon displacements. 
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FIG. 2. (Color online) Bias dependence of the phonon- 
mediated conductance G p h- (a) G p h for different tempera- 
tures T, _Ef=0.26 eV and = 21.8°. From top to bottom, 
T = 500 K (light green), 400 K (green), 300 K (red), 200 K 
(blue), 100 K (magenta), 60 K (cyan), 20K (black). The black 
dashed curves are single parameter fits described in text. This 
scaling breaks down at temperatures below the energy of the 
beating phonon mode, (b) G v h at T = 20 K and Ef=Q-2Q eV 
for three different structures: from top to bottom, 9 = 13.2° 
(green curve), 9 = 21.8° (red curve), and 9 = 27.8° (blue 
curve). The width of the plateau reflects the beating phonon 
energy (Fig. [3b). 

The results for the phonon mediated conductance arc 
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summarized in Fig. [T] and compared with the electronic 
conductances and band splittings. We include all the 
commensurate structures up to N atoms = 4L 2 /a 2 ps 1100 
atoms per rotational supercell. The splitting A is ob- 
tained from the low-energy effective 4x4 Hamiltonian 
[3 [lj| calculated using Eq. ©. The specific commen- 
surate structure is_characterized by a graphene lattice 
IS" 



vector 17 



21 



so a priori it need not follow a 
smooth curve when plotted against any single scalar vari- 
able. In Fig. [2(a) and (b) the splitting appears chaotic 
when plotted against angle. However, a simple smooth 
trend emerges when plotted against the size of the super- 
cell, reflecting the role of umklapp processes 
As expected, in Fig QJc) and (d) we see generally sim- 
ilar behavior for the calculated electronic conductance 

Gel 
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Both A and G e i decay with increasing supercell 
size, such that G e i decreases by 7 orders of magnitude 
when L increases from 2 nm to 6 nm, which is typi- 
cal of the range seen in bilayer graphene Moire patterns 
23l 35l - 39J . If phonon scattering is omitted, the calcula- 
tion of electronic conductance requires some phenomeno- 
logical broadening, and is proportional to the value of this 
broadening [l9| . We therefore give the electronic conduc- 
tance relative to that of the smallest nontrivial rotated 
structure, i.e. since the ratio is independent 

of the broadening. G7 was estimated to be G-j w 30 
mS//im 2 for E F = 0.26 eV. 

When phonon scattering is included, the behavior is 
radically different in two distinct ways. First, there is a 
smooth dependence on rotation angle and not on super- 
cell size, the opposite of the behavior without phonons. 
Second, the maximum variation with angle (excluding 
ideal AB stacking and very small angles) is only about a 
factor of 10, rather than the exponential variation seen 
without phonons in Fig. Q}i. Thus phonon-mediated con- 
duction is expected to dominate except perhaps near ori- 
entations where the supercell is particularly small. 

The exponential decay of electronic coupling A with 
supercell size L, shown in Fig. [T]b, can be understood 
from previous work firMTol [25I ] . As discussed there, states 
near the K point of one layer are coupled to the rotated K 
point of the other layer by Umklapp process. The larger 
the commensuration period L, the larger the reciprocal 
lattice vector Q needed for Umklapp coupling. The cou- 
pling is related to the Q Fourier component of the finite- 
range interlayer hopping amplitude. The conductance in 
turn scales as A 2 , Fig. [TJi. The solid line in Fig. QJi is 
G e i oc exp(— L/lei) with the decay length l e i ~ 0.25 A. 
However, a few points show striking deviations, by up to 
a factor of 50. This is because the electronic contribution 
is very sensitive to translation d [l9[ as well as rotation 
24\. 

In contrast to a single graphene layer, the bilayer has 
low-energy optical flexural phonons that can easily pro- 
vide the momentum needed to couple the rotated K 
points without Umklapp, so there is no exponential de- 



pendence on periodicity in Fig.QJf) — indeed, no obvious 
dependence on periodicity at all. The full numerically 
results of Fig. [lje) can be understood using a simplified 
model |24j which considers coupling between the rotated 
Dirac points by the beating-mode phonon: 



G ph oc g E F 



2 n 9 



UJ() 



(4) 



where ujg = LOb{<lK) is the frequency of the beating mode 
ujb at wavevector qx connecting the rotated Dirac points, 
ng is the thermal population of ujg phonon, and the con- 
stant of proportionality is given in Ref. [24j. This model 
has only one adjustable parameter, the effective electron- 
phonon matrix element g. As expected, in Fig. QJe) the 
simple model breaks down for small rotation angles, be- 
cause can no lon ger be neglected relative to qx- 
Previous work lj|, 25 1 emphasized that the purely elec- 
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FIG. 3. (Color online) (a) Relative contribution of phonons 
at different energies to the conductance. Peaks in G u clearly 
indicate the dominant role of specific phonons. (b) Energy of 
beating phonon mode at the wavevector qx connecting the 
rotated Dirac points is shown as a function of rotation angle 
9. The dashed curve is a fit described in the text. As the 
angle increases, the relevant phonon energy increases. The 
dashed lines show that for each angle, this phonon energy 
corresponds well to the structure in (a) for that angle. Inset 
schematically illustrates beating mode geometry, (c) Tem- 
perature dependence of the low bias conductance in the same 
structures. 
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tronic conductance shows an extraordinary sensitivity to 
fine details of the model used in calculating the matrix 
elements. This occurs because subtle changes in the an- 
alytic form of the interlayer interactions can radically 
change the asymptotic behavior at large Q, which con- 
trols Umklapp scattering. Once we include phonon scat- 
tering, this extraordinary sensitivity disappears, as does 
the strong sensitivity to lateral translation. (Details and 
further calculations are given in Ref. 0.) Most impor- 
tantly, phonons not only solve the "numerical" problems 
of hypersensitivity, but also solve the conceptual prob- 
lems that arise in the small-scattering limit [19J , so we can 
calculate absolute conductance without any phenomeno- 
logical broadening. 

The conductance depends sensitively on temperature, 
and also depends on bias voltage, as shown in Fig. [2^. 
The behavior is particularly interesting at low temper- 
ature, Fig. [2Jd. At sufficiently low temperatures, the 
phonon-mediated current shows a diode-like behavior, 
turning on at a sharply defined threshold voltage. The 
threshold voltage Vth depends on the rotation angle. 
The sharpness of the current onset depends on the ra- 
tio k B T/eV t h- 

The general behavior at high biases can be understood 
from simple density of states argument (see Eq. (jSllO . 
which suggests scaling in the limit Ep 3> ksT [281 ] : 
G ph oc (eV)' 1 J^_ eV E 2 dE = E% - eVE F + e 2 V 2 /3. 



As shown in Fig. [2jt this scaling works extremely well 
at high temperatures (dashed curves), reproducing the 
full numerical result with only one adjustable parame- 
ter chosen to match the zero-bias conductance at each 
temperature. 

The current is suppressed at low temperatures and low 
biases, because the required phonons cannot be excited 
either thermally or electronically. As the bias becomes 
large enough for the carriers to excite a phonon, conduc- 
tance starts to increase as shown in Fig. (A some- 
what related effect has been suggested in tunneling from 
a metal tip to graphene 4(], 41 1.) The relevant phonon 
energy uig increases with the rotation angle, Fig. [3}d, so 
a larger threshold bias is needed to allow phonon as- 
sisted tunneling as seen in Fig. Ox As a result, the 
low-conductance plateau gets wider as Vth oc ojg and the 
"on/off" conductance ratio increases as oc 1/ne- 

To see clearly which phonons contribute most to the 
total conductance, in Fig. [3ji we plot the fraction of con- 
ductance arising from phonons at a specific energy. In 
all three cases, the energy of the main peak corresponds 
to the beating mode connecting the rotated Fermi cir- 
cles. This is shown by the dashed lines in Fig. [5b, those 
phonon energies at 95 cm -1 , 150 cm -1 , 220 cm -1 roughly 
coincide with the onsets of the phonon contribution to the 
conductance for the three structures calculated in Fig.[3ji. 
(The in-phase flexural phonon, whose energy vanishes at 
the r-point, does not contribute to the conductance.) 

The dashed line in Fig. [3Jd shows a fit to the beat- 



ing phonon mode phenomenological dispersion, u>b{<lK) = 

1 /2 

\oj 2 + Qk k /p\ > where ur is zone center beating 
phonon frequency, p is the graphene mass density, and 
the single fitting parameter is the effective bending stiff- 
ness k [42I ]. The energy wt, (<?#-) of the phonon cou- 
pling Dirac points also sets the temperature scale for the 
crossover from a linear to a sublinear behavior, as shown 
in Fig. [3b. 

At low rotation angles, the phonon-mediated conduc- 
tance is expected to be larger than that in the calcula- 
tions above due to the electronic coupling to the inter- 
layer slidin g m ode [43| , and to a ny Fermi velocity renor- 
malization [20j| . We estimate [24j that the current is 
underestimated by less than 10% for angles larger than 
9 > 10°. 

In conclusion, we have shown that phonon scattering is 
a dominant factor in the interlayer conduction of twisted 
bilayers of graphene. It relaxes the momentum conser- 
vation, eliminating the extraordinary sensitivity to ro- 
tation and translation of the layers, and it drastically 
changes the functional dependence. Even at low temper- 
ature, phonons can dominate the conductance, provided 
that the bias across the layers is large enough to allow 
the excitation of a phonon during the tunneling. This 
effect leads to intriguing diode-like electrical characteris- 
tics with high "on/off" conductance ratio, provided that 
other scattering mechanisms are small. These phenom- 
ena should be relevant in other two dimensional inter- 
faces where current is suppressed by lateral momentum 
mismatch of the electronic states, which may be expected 
in novel tunneling devices 44 1 . 
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ANALYTICAL EXPRESSION FOR G ph 

The full numerically results of Fig. le of the main text 
can be understood using a simplified model. The phoiion 
must provide the momentum connecting the respective 
Fermi circles, q K « 2Ks'm(9/2). Here K = 47r/(3a) ps 17 
nm _1 , and we neglect the Fermi circle radius kp rel- 
ative to the separation qx- The dominant contribu- 
tion comes from the beating mode (i.e. the out-of- 
phase bilayer flexural mode). Approximating the asso- 
ciated electron-phonon matrix element as momentum- 
independent, \(n\H^ ph \H)\ 2 « g 2 h/ (2M c u; g ), where 
tug = ojb(qj() is the frequency of the beating mode u>b at 
wavevector qK- Then the low bias conductance can be 
evaluated from Eqs. (1) and (2) of the main text: 

G ph = 2TTe 2 g 2 neE F A /{McuJ8TT 2 fr 4 v F ) (SI) 

where Mq is the carbon mass, Aq = a 2 \/3/2 is the 
area of the two atom cell, vf = y/3ta/(2h) ~ 10 s cm/s 
is the Fermi velocity, ng is the thermal population of 
cog phonon, and g is an effective (momentum averaged) 
electron-phonon constant. With a single adjustable pa- 
rameter g ps 0.34 eV/A, Eq. (fSTj) accounts very well for 
the numerical result over a wide range of rotational an- 
gles as shown by the solid curve in Fig. 1(e) of the main 
text. As expected, the simple model breaks down for 
small rotation angles, because \zf can no longer be ne- 
glected relative to qx- 



At low rotation angles, the phonon-mcdiatcd con- 
ductance is expected to be larger than that in the 
calculations above due to the electronic coupling to 
the interlayer sliding mode, and to any Fermi velocity 
renormalization[20j . First, we estimate the sliding mode 
contribution by assuming the electron-phonon coupling 
strength for the sliding mode and beating phonons are 
similar ~ . Using acoustic sound velocity of v s — 

20 x 10 3 m/s and the zone center sliding mode energy in 
the bilayer [43| hu> s o = 31 cm -1 we can estimate the 
relevant phonon frequency as w s = \A^ + v 2 s q K . Ac- 
cording to Eq. (|Sip . the total conductance would scale as 
G p h(l +ujg/uj 2 ). This increases the conductance by only 
7% for 6 = 10°, increasing to 25% and 65% for 6 = 5° and 
3° respectively. At the doping level used in Fig. 1 of the 
main text, the respective Fermi circles start to overlap 
at 6 < 3°. To estimate the Fermi velocity renormaliza- 
tion v F , we use Ref. (20|: vp/v F = 1 — 9[t±/(HvpqK)] 2 , 
where t± = 0.11 cV. The Fermi energy Ep = hvpy/nn s 
used in Eq. (|S1[) will also be reduced for a fixed concen- 
tration n s = 5 x 10 12 cm~ 2 , increasing the conductances 
by roughly a factor of vp/vp ps 2.14, 1.26, and 1.06 for 
(9 = 3° , 5°, and 10° respectively. Therefore, our results 
in Fig. le and If of the main text are expected to be 
quantitatively robust (within 10%) for angles larger than 
6 > 10°. At smaller angles the phonon-mediated conduc- 
tance is larger than calculated here, further enhancing 
the interlayer transport. 



ADDITIONAL CONTRIBUTIONS TO G ph AT 
LOW ROTATION ANGLES 
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FIG. 4. (Color online) LDA calculations of the band splitting 
A in AB stacked bilayer graphene. A versus the interlayer 
distance h: from LDA calculations (blue squares), from the 
model Hamiltonian Eq. (3) in the main text (red hollow cir- 
cles), and an exponential fit: A = Aoexp((/io — h)/\) with 
A = 0.6 A and A = 0.35 eV at ho = 3.35 A. 
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FIG. 5. (Color online) Tight-binding model test against LDA 
results for the velocity renormalization. The blue squares are 
results from Ref. [2l| in structures from left to right: (n, m) = 
(6,7), (4, 5), (5, 9), (1, 3), and (1, 4). The hollow red circles 
are our tight-binding results in structures from left to right: 
(19, 20), (15, 16), (12, 13), (10,11), (8, 9), (6, 7), (4, 5), (5, 
9), (1, 3), (1, 4), (1, 10), (1, 20), (1, 35), (1, 55). The black 
solid curve gives a perturbation theory result from Ref. [20T ] : 



Vf/vf 



)[t±/(hv F qK)] , where t± = 0.11 eV. 
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FITTING OF OUR TIGHT-BINDING MODEL 
TO LDA RESULTS 



To obtain reasonable quantitative accuracy, we de- 
sire a model that correctly describes the dependence on 
both rotation and interlayer spacing, insofar as these 
are known from LDA calculations. We require that the 
model reproduce the LDA results for the dependence of 
A on layer spacing, which we calculated in the Bernal 
stacked bilayer as shown in Fig. @] Our model in Eq. (3) 
of the main text contains four parameters: t±, X Zl X xy 
and a. We have used only three LDA values to fit them: 
A in Bernal stacked bilayer, interlayer spacing depen- 
dence in Bernal stacked bilayer (i.e. derivative dA/dh at 
h = Hq), and A value in 7 x 7 twisted bilayer. The model 
could be further improved by fitting to more LDA calcu- 
lations. Nevertheless, we tested that our model describes 
reasonably well other LDA calculations, such as Fermi 
velocity renormalization shown in Fig. [5j which were not 



used in fitting the parameters of the tight-binding model. 

We emphasize that the qualitative points of our paper 
could have been obtained with the previous models. We 
introduce a new model simply as a tool to get reason- 
able accuracy for the electron-phonon coupling, since the 
previous models were not designed to give the correct in- 
terlayer distance dependence. In Figs. [6] and [7] we show 
that variation of the parameters of the model by as much 
as a factor of two does not change the phonon mediated 
conductance by more than 30%. This can be contrasted 
with the electronic conductance, which is far much more 
sensitive to the parameters of the tight-binding model. 
The same is true of the sensitivity to lateral sliding of 
one layer with respect to the other, also shown in Figs. [6] 
and [71 

Finally, in Fig. [5] shows that the long range nature of 
the tight-binding model in Eq. (3) of the main text in- 
troduces relatively small electron-hole asymmetry, so all 
the results in the main text are reported for the electron 
branch. 
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FIG. 6. (Color online) Sensitivity of electronic conductance 
G e i and phonon mediated conductance G p h to the lateral cut- 
off distance X xy and on the sliding shift d. Here Ao = 1.7 
A and the sliding shift d (in units of carbon bondlength 
Rcc) is along the C-C bond direction in the AA stacked bi- 
layer after rotation, with 6 = 21.8°. (a) and (b) show G e i at 
T — K, while (c) and (d) show G p h at room temperature 
and Ef = 0.26 eV. Note, that d = corresponds to the struc- 
ture where AA stacked bilayer is rotated around AA axes and 
translational shift d — 3Rcc gives an equivalent structure. 
Conventionally, twisted structures are obtained from the AB 
stacked bilayer rotated around the AB axes by a given angle, 
which would correspond to structure with d = 2Rcc used 
for calculations in the main text. Structure with d — Rcc 
corresponds to the conventional structure rotated by a com- 
plimentary angle 60° — 6, i.e. 6 = 38.2° 
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FIG. 7. (Color online) Sensitivity of electronic conductance 
Gel and phonon mediated conductance G p h to the value of A. 
Here Ao = 0.6 A, see Eq. (3) in the main text. (a) G e i at 
T = K. (b) Gph at room temperature and Ef = 0.26 eV. 
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FIG. 8. (Color online) Fermi energy dependence of the 
phonon mediated conductance. G p h in 9 fts 21.8° bilayer 
structure at T=300 K (red circles). The black solid curve is a 
parabolic fit motivated by the approximation of momentum 
independent electron-phonon matrix elements, i.e. Eq. (|S1I) . 
The single parameter fit does not capture the slight electron- 
hole asymmetry which comes from the long-range nature of 
the interaction in Eq. (3) of the main text. The small resid- 
ual conductance at Ef = is due to the thermal electron-hole 
pair excitations. 



